Application: Mincer Wage Equation

Mincer wage equation (Mincer 1958):

\[\ln(wage)_i = \beta_0 + \beta_1 \mbox{exp}_i + \beta_2 \mbox{exp}_i^2 + e_i\]

• For ease of exposition, dropping education.

Application: Mincer Wage Equation

Mincer wage equation (Mincer 1958):

\[\ln(wage)_i = \beta_0 + \beta_1 \mbox{exp}_i + \beta_2 \mbox{exp}_i^2 + e_i\]

• Average effect of experience: \[\beta_1 + 2 \cdot \beta_2 \cdot \mathbb{E}(\mbox{exp}),\]

• Experience level that maximizes log wage: \[\left | \frac{\beta_1}{2 \beta_2}\right |.\]

Application: Mincer Wage Equation

• Let \(\widehat \beta_{1,N}, \widehat \beta_{2,N}\) denote OLS estimators,

• Let \(\overline{\mbox{exp}}_N\) denote sample mean of experience.

• Estimate average effect of experience by \[\widehat \beta_{1,N} + 2 \cdot \widehat \beta_{2,N} \cdot \overline{\mbox{exp}}_N,\]

  • Estimator is a non-linear function of \((\widehat \beta_{1,N}, \widehat \beta_{2,N}, \overline{\mbox{exp}}_N)\).

Application: Mincer Wage Equation

• Let \(\widehat \beta_{1,N}, \widehat \beta_{2,N}\) denote OLS estimators,

• Estimate experience level that maximizes log wage by \[\left | \frac{\widehat \beta_{1,N}}{2 \widehat \beta_{2,N}}\right|.\]

  • Estimator is a non-linear functions of \((\widehat \beta_{1,N}, \widehat \beta_{2,N})\).

Application Mincer Wage Equation

library (AER)
data(CPS1985)
df  <- CPS1985[ CPS1985$education==12, c("wage","experience")]
mean.exp <-mean(df$experience) 
mean.exp
reg.1 <- lm(I(log(wage))~experience+I(experience^2),data=df)
reg.1
a <- c(0,1,2*mean.exp)
avg.eff.exp <- as.numeric(t(a)%*%coef(reg.1)) 
maxexp <- abs(coef(reg.1)[2]/(2*coef(reg.1)[3]))

• Consider estimation using CPS 1985 data.

  • available in package AER.

  • after loading AER, see ?CPS1985 for description of data.

• Consider only using observations with \(12\) years of schooling.

Application Mincer Wage Equation

library (AER)
data(CPS1985)
df  <- CPS1985[ CPS1985$education==12, c("wage","experience")]
mean.exp <-mean(df$experience) 
mean.exp
reg.1 <- lm(I(log(wage))~experience+I(experience^2),data=df)
reg.1
a <- c(0,1,2*mean.exp)
avg.eff.exp <- as.numeric(t(a)%*%coef(reg.1)) 
maxexp <- abs(coef(reg.1)[2]/(2*coef(reg.1)[3]))

[1] "Mean Experience: 18.14"

Call:
lm(formula = I(log(wage)) ~ experience + I(experience^2), data = df)

Coefficients:
    (Intercept)       experience  I(experience^2)  
      1.5829802        0.0352680       -0.0005831  

Application Mincer Wage Equation

mean.exp <-mean(df$experience) 
mean.exp
reg.1 <- lm(I(log(wage))~experience+I(experience^2),data=df)
reg.1
a <- c(0,1,2*mean.exp)
avg.eff.exp <- as.numeric(t(a)%*%coef(reg.1)) 
maxexp <- abs(coef(reg.1)[2]/(2*coef(reg.1)[3]))
avg.eff.exp
maxexp

[1] "Estimated Avg Effect Experience: 0.0141"

[1] "Estimated Experience that Maximizes Wage: 30.2"

Application Mincer Wage Equation

mean.exp <-mean(df$experience) 
mean.exp
reg.1 <- lm(I(log(wage))~experience+I(experience^2),data=df)
reg.1
a <- c(0,1,2*mean.exp)
avg.eff.exp <- as.numeric(t(a)%*%coef(reg.1)) 
maxexp <- abs(coef(reg.1)[2]/(2*coef(reg.1)[3]))
avg.eff.exp
maxexp

Application: Mincer Wage Equation

• How to compute standard errors? construct confidence intervals? perform hypothesis tests?

Traditional approach: Delta Method.

• Asymptotic approximation based on linear approximation from Mean Value Theorem combined with LLN, CLT.

• Limitations:

  • only justified asymptotically,
  • sometimes provides poor approximation,
  • derivations using delta method can be tedious, chance of algebraic errors.

Bootstrap principle

• Treat sample of size \(N\) as if it were the population.

• Perform MC simulation, creating many bootstrap samples of size \(N\) by resampling with replacement from original sample.

• Compute estimator on each sample, use distribution of estimates across samples to approximate distribution of estimator.

• Difference between drawing from population vs from sample should vanish as \(N \rightarrow \infty\).

Bootstrap: Basic Principle

Bootstrap Method

• Limitations:

  • only justified asymptotically,
  • sometimes provide poor approximation,
  • can be computationally intensive

• Advantages:

  • Can bypass derivations, complicated formulas. Less tedious to implement. Removes chance of algebraic error in derivations.
  • Sometimes provides better approximation than delta method.

Bootstrap Method

• There are many versions of the bootstrap, but we will focus on nonparametric bootstrap and corresponding . . .

  • estimate of bias,
  • standard errors,
  • normal-approximation CI,
  • percentile CI,
  • percentile-t CI,
  • hypothesis tests.

Sample and Estimator

• Suppose \(\mathbf{W} = (W_1,W_2,...,W_N)\) is observed sample
  • e.g., \(W_i = (Y_i,X_i)\) in regression context.
• Estimator \(\widehat \theta = s(\mathbf{W})\)
  • expressing estimator as function of data, rule for how to use data to compute estimate.
  • in our application, \(\widehat \theta\) is estimator of average effect of experience or of experience level that maximizes wages.

Example: Bootstrap Estimate Bias

• Given \(B\) bootstrap replications, define the boostrap estimate of bias by \[\widehat{\mbox{bias}}^*= \frac{1}{B} \sum_{b=1}^B \widehat \theta^*(b) - \widehat \theta \] where \(\widehat \theta\) is estimate from original sample.

• why centered at \(\widehat \theta\)?

Example: Algorithm for Percentile CI

• Given \(B\) bootstrap replications, define \(q^*_{\alpha}\) as the \(\alpha\) empirical quantile of \((\widehat \theta^*(1),\widehat \theta^*(2),…,\widehat\theta^*(B))\)

  • compute \(q^*_{\alpha}\) as the value such that fraction \(\alpha\) of the bootstrap replications are less than \(q^*_{\alpha}\).

• Bootstrap Percentile CI: \[C^{PC}= [q_{\alpha/2},~ q_{1-\alpha/2}].\]

Example: Algorithm for Percentile-t CI

• Given \(B\) bootstrap replications, define \(q^*_{\alpha}\) as the \(\alpha\) empirical quantile of \((T^*(1), T^*(2), ...,T^*(B))\), where \[T^*(b) = \frac{\widehat \theta^*(b) - \widehat \theta}{s(\widehat \theta^*(b))}\]

where \(s(\widehat \theta^*(b))\) is the standard error of \(\widehat \theta^*(b)\) computed using only the \(b\)th bootstrap replication.

• Note centering at \(\widehat \theta\), and note s.e. using only \(b\)th bootstrap replication.

Example: Algorithm for Percentile-t CI (cont’d)

• Given \(B\) bootstrap replications, define \(q^*_{\alpha}\) as the \(\alpha\) empirical quantile of \((T^*(1), T^*(2), ...,T^*(B))\).

• Bootstrap Percentile-t CI: \[C^{Pt}= [\widehat \theta - s(\widehat \theta) q_{1-\alpha/2},~ \widehat \theta - s(\widehat \theta) q_{\alpha/2}].\]

Pros and Cons Bootstrap Methods

• Bootstrap standard errors:
  • Pros:
    • intuitive, easy to implement,
    • researchers familiar with s.e., is the norm to report s.e..
  • Cons:
    • Valid only under relatively strong conditions,
    • May be erratic, poorly behaved.

Pros and Cons Bootstrap Methods

• Normal-approximation bootstrap CI:

  • Pros:
    • Intuitive, easy to implement,
    • Researchers familiar with normal CIs,
    • Many estimators are asymptotically normal.

Pros and Cons Bootstrap Methods

• Normal-approximation bootstrap CI:

  • Cons:
    • Invalid if bootstrap s.e. are invalid,
    • Poorly behaved if bootstrap s.e. are poorly behaved,
    • Normality may be bad approximation
      • especially poor approximation if estimator is biased and/or has highly asymmetric distribution.

Pros and Cons Bootstrap Methods

• bootstrap percentile CI:

  • Pros:

    • intuitive, easy to implement,

    • valid under weaker conditions than bootstrap s.e. or bootstrap normal CI.

  • Cons:

    • can provide poor approximation if estimator is biased and/or has highly asymmetric distribution.

Pros and Cons Bootstrap Methods

• bootstrap percentile-t CI:

  • Pros:

    • valid under weaker conditions than bootstrap s.e. or bootstrap normal CI.

    • better theoretical properties (more accurate) than percentile CI.

  • Cons:

    • need to compute standard error on each bootstrap replication.

Avg Effect of Experience

• We now implement bootstrap s.e., CIs and hypothesis tests for average effect of experience.

  • \(\theta = \beta_{1} + 2 \cdot \beta_{2} \cdot \mathbb{E}[\mbox{exp}],\)

  • \(\widehat \theta = \widehat \beta_{1,N} + 2 \cdot \widehat \beta_{2,N} \cdot \overline{\mbox{exp}}_N.\)

• We first proceed in base R, then consider using package boot.

• First step: create bootstrap samples.

Avg Effect of Experience

• Create bootstrap sample:

  • Randomly draw \(N\) integers, with replacement, from \((1,...,N)\), with each integer equally likely to be draw,

  • Create new sample formed by observations with those indices.

  • Some observations from original sample will be drawn multiple times into bootstrap sample, others will not appear in bootstrap sample.

N <- nrow(df) 
indices.b1<-sample(1:N,replace=TRUE)  # draw indices with replacement from (1,...,N), N sample size
tab.b1 <- df[indices.b1,] # construct bootstrap sample

Bootstrap Replications Avg. Effect Experience

• given bootstrap sample, \(\mathbf{W}^{*b}=(W^{*b}_1,W^{*b}_2,…,W^{*b}_N)\), compute \(\widehat \theta^*(b)=s(\mathbf{W}^{*b}).\)

  • compute mean experience on bootstrap sample, \(\overline{\mbox{exp}}^{~*b}\).
  • run regression on bootstrap sample, obtain bootstrap sample OLS coefficients \((\widehat \beta_{0}^{*b},\widehat \beta_{1}^{*b}, \widehat \beta_{2}^{*b})\)
  • compute \[\widehat \theta^*(b) = \widehat \beta_{1}^{*b} + 2 \cdot \widehat \beta_{2}^{*b} \cdot \overline{\mbox{exp}}^{~*b}.\]

Bootstrap Replications Avg. Effect Experience

#  tab.b1 <- df[indices.b1,], bootstrap sample
mean.exp.b1 <- mean(tab.b1$experience) #  sample mean experience on bootstrap sample
reg.b1 <- lm(I(log(wage))~ experience+I(experience^2),data=tab.b1) # OLS on bootstrap sample
# form bootstrap replication of estimated mean effect experience
a <- c(0,1,2*mean.exp.b1)
theta.b1 <- as.numeric(t(a)%*%coef(reg.b1)) 

Construct \(B\) Bootstrap Replications

• We have created 4 bootstrap samples, and corresponding 4 bootstrap replications.

• We want to create \(B\) bootstrap samples, and corresponding \(B\) bootstrap replications, for \(B\) large, for example, \(B=10,000\).

• How to repeat steps many times?

  • we first consider approach in base R using replicate,
  • we then consider using functions from package boot.

Use replicate to Construct \(B\) Boot. Replications

f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
N <- nrow(df)  # Number of obs in sample
B <- 10000 # number of bootstrap replications
boot.est <-replicate(B, f.boot1(df, sample(1:N, N, replace=TRUE)))

• Create user defined function and use replicate to generate \(B\) bootstrap samples, and corresponding \(B\) bootstrap replications.

Use Replicate to Construct \(B\) Boot. Replications

f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
N <- nrow(df)  # Number of obs in sample
B <- 10000 # number of bootstrap replications
boot.est <-replicate(B, f.boot1(df, sample(1:N, N, replace=TRUE)))

• User defined function takes two arguments:
  1. original sample as data frame,
  2. indices;
• Function returns bootstrap replication.

Use Replicate to Construct \(B\) Boot. Replications

f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
N <- nrow(df)  # Number of obs in sample
B <- 10000 # number of bootstrap replications
boot.est <-replicate(B, f.boot1(df, sample(1:N, N, replace=TRUE)))

• Inputting user defined function, original data and random sample of indices into replicate function, which is returning \(B\) bootstrap replications.

• Use output to compute bootstrap s.e. or CI.

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1

• user defined function takes two inputs:
  • data frame (original sample)
  • indices;
• function returns bootstrap replication.

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1

• boot taking as inputs:

  • original sample (data frame),
  • user defined function,
  • number of replications.

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1

• to understand output from boot, enter:

  • ?boot
  • str(boot.out1)

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = df, statistic = f.boot1, R = B, stype = "i")


Bootstrap Statistics :
      original        bias    std. error
t1* 0.01410963 0.00004975603 0.002904684

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1

• boot.ci taking as inputs:

  • saved output from boot,
  • choice of confidence level,
  • type of bootstrap CI.

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1

• to understand output from boot.ci, enter:

  • ?boot.ci
  • str(boot.ci1)

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.out1, conf = 0.95, type = c("norm", "perc", 
    "bca"))

Intervals : 
Level      Normal             Percentile            BCa          
95%   ( 0.0084,  0.0198 )   ( 0.0085,  0.0199 )   ( 0.0083,  0.0197 )  
Calculations and Intervals on Original Scale

Analytical s.e.

Let \(\sigma^2_X = \mbox{Var}(X_i)\), and

\[ a = ~ \left[ \begin{array}{c} 0\\ 1 \\ 2 \mathbb{E}(\mbox{exp}) \end{array} \right]^{\prime}, \qquad\qquad \widehat a= \left[ \begin{array}{c} 0\\ 1 \\ 2 \overline{\mbox{exp}} \end{array} \right]^{\prime}. \]

Analytical s.e. ignoring sampling variability in \(\overline{\mbox{exp}}\), \[\mbox{s.e.}(\widehat \theta) = \frac{1}{\sqrt{N}} \left( \widehat a^{\prime} \widehat \Sigma \widehat a \right)^{1/2}\]

Analytical s.e.

Let \(\sigma^2_X = \mbox{Var}(X_i)\), and \[ a = ~ \left[ \begin{array}{c} 0\\ 1 \\ 2 \mathbb{E}(\mbox{exp}) \end{array} \right]^{\prime}, \qquad\qquad \widehat a= \left[ \begin{array}{c} 0\\ 1 \\ 2 \overline{\mbox{exp}} \end{array} \right]^{\prime}. \]

Analytical s.e. using delta method, \[\mbox{s.e.}(\widehat \theta) = \frac{1}{\sqrt{N}} \left( \widehat a^{\prime} \widehat \Sigma \widehat a + 4 \widehat \beta_2^2 \widehat \sigma^2_X \right)^{1/2}.\]

Analytical s.e.

# not accounting for variability in sample mean
library(sandwich)
a <- c(0,1,2*mean.exp)
se.1 <-  as.numeric(sqrt(t(a)%*%vcovHC(reg.1)%*%a))
signif(se.1,3)

[1] 0.00277

# corresponding 95% CI
c(signif(theta.hat -qnorm(.975)*se.1,3), signif(theta.hat+qnorm(.975)*se.1,3))

[1] 0.00869 0.01950

# using delta method
se.2 <- as.numeric(sqrt(t(a)%*%vcovHC(reg.1)%*%a +
                          4*(coef(reg.1)[[3]])^2*var(df$experience)/nrow(df)))
signif(se.2,3)

[1] 0.00293

# corresponding 95% CI
c(signif(theta.hat -qnorm(.975)*se.2,3), signif(theta.hat +qnorm(.975)*se.2,3))

[1] 0.00836 0.01990

Bootstrap Percentile-t CI

• For percentile-t CI, need to estimate variance of estimator on each bootstrap sample.

• Can do so using result from Delta method.

Bootstrap Percentile-t CI

f.boot1.t <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  est.b <- as.numeric(t(a)%*%coef(reg.b)) 
  # form estimated variance
  var.b <- as.numeric(t(a)%*%vcovHC(reg.b)%*%a+4*(coef(reg.b)[[3]])^2*var(d[i,]$experience)/nrow(d[i,]) )
  c(est.b,var.b)
  }
boot.out1.t<- boot(data=df, statistic=f.boot1.t, R=B, stype="i")

• Notice that function f.boot1.t is returning variance of estimator computed on bootstrap sample as second component, necessary if inputting into boot for later input into boot.ci to construct percentile-t CI.

Bootstrap Percentile-t CI

f.boot1.t <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  est.b <- as.numeric(t(a)%*%coef(reg.b)) 
  # form estimated variance
  var.b <- as.numeric(t(a)%*%vcovHC(reg.b)%*%a+4*(coef(reg.b)[[3]])^2*var(d[i,]$experience)/nrow(d[i,]) )
  c(est.b,var.b)
  }
boot.out1.t<- boot(data=df, statistic=f.boot1.t, R=B, stype="i")
summary(boot.out1.t$t)

       V1                 V2             
 Min.   :0.002754   Min.   :0.000004574  
 1st Qu.:0.012117   1st Qu.:0.000007711  
 Median :0.014067   Median :0.000008607  
 Mean   :0.014086   Mean   :0.000008716  
 3rd Qu.:0.016044   3rd Qu.:0.000009618  
 Max.   :0.025135   Max.   :0.000015271  

V2 is estimated variance on bootstrap sample.

Bootstrap Percentile-t CI

boot.out1.t<- boot(data=df, statistic=f.boot1.t, R=B, stype="i")
boot.ci1.t <- boot.ci(boot.out1.t, conf=0.95, type=c("norm","perc","bca","stud"))
boot.ci1.t

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.out1.t, conf = 0.95, type = c("norm", 
    "perc", "bca", "stud"))

Intervals : 
Level      Normal            Studentized     
95%   ( 0.0083,  0.0199 )   ( 0.0084,  0.0199 )  

Level     Percentile            BCa          
95%   ( 0.0085,  0.0199 )   ( 0.0084,  0.0199 )  
Calculations and Intervals on Original Scale

Rerun Analysis?

• Can rerun bootstrap simulation to see if results are stable across simulations.
• If not?
  • insufficient number of replications? increase \(B\)?
  • problem with bootstrap method in particular example?
    • sometimes bootstrap standard errors invalid/unstable, arises when bootstrap replications take extreme values, e.g., when dividing by number close to zero for some replications.

Rerun Analysis?

• Can rerun bootstrap simulation to see if results are stable across simulations.
• If not?
  • problem with bootstrap method in particular example?
    • sometimes bootstrap standard errors invalid/unstable, arises when bootstrap replications take extreme values, e.g., when dividing by number close to zero for some replications.
    • percentile CIs are valid under weaker assumptions and tend to be more robust/stable.

Implement: Exp Level Max. Wages

• We now implement bootstrap s.e. and CIs for experience level that maximizes log wage.

  • \(\theta = \left | \frac{\beta_{1}}{2 \cdot \beta_2} \right|,\)
  • \(\widehat \theta = \left | \frac{\widehat \beta_{1}}{2 \cdot \widehat \beta_2} \right|,\)

• Key difference from before:

  • parameter involves ratio of coefficients,
  • estimator divides by estimated coefficient,
  • cause for worry?

Max Experience

f.boot2 <- function(d,i){ 
   reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
   abs(coef(reg.b)[2]/(2*coef(reg.b)[3])) }
boot.out.b<- boot(data=df, statistic=f.boot2, R=B, stype="i")
boot.out.b
boot.ci.b <- boot.ci(boot.out.b, conf=0.95, type=c("norm","perc","bca"))
boot.ci.b

Max Experience

f.boot2 <- function(d,i){ 
   reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
   abs(coef(reg.b)[2]/(2*coef(reg.b)[3])) }
boot.out.b<- boot(data=df, statistic=f.boot2, R=B, stype="i")
boot.out.b
boot.ci.b <- boot.ci(boot.out.b, conf=0.95, type=c("norm","perc","bca"))
boot.ci.b

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = df, statistic = f.boot2, R = B, stype = "i")


Bootstrap Statistics :
    original   bias    std. error
t1* 30.23942 4.915459    217.3218

Max Experience

f.boot2 <- function(d,i){ 
   reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
   abs(coef(reg.b)[2]/(2*coef(reg.b)[3])) }
boot.out.b<- boot(data=df, statistic=f.boot2, R=B, stype="i")
boot.out.b
boot.ci.b <- boot.ci(boot.out.b, conf=0.95, type=c("norm","perc","bca"))
boot.ci.b

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.out.b, conf = 0.95, type = c("norm", 
    "perc", "bca"))

Intervals : 
Level      Normal             Percentile            BCa          
95%   (-449.04,  498.47 )   (  25.14,   52.12 )   (  25.22,   53.14 )  
Calculations and Intervals on Original Scale

Bootstrap Percentile-t CI

f.boot2.t <- function(d,i){ 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  est.b <- abs(coef(reg.b)[2]/(2*coef(reg.b)[3]))  
  a <- c(0,-1/(2*coef(reg.1)[[3]]),coef(reg.1)[[2]] /(2*(coef(reg.1)[[3]])^2))
  var.b <- as.numeric(t(a)%*%vcovHC(reg.1)%*%a)
  c(est.b,var.b)
  }
boot.out2.t<- boot(data=df, statistic=f.boot2.t, R=B, stype="i")

Test Using Asymptotic Normality

Consider

• \(H_0: ~ \theta = \theta_0,\)
• \(H_1: ~ \theta \ne \theta_0\).

Test based on asymptotic normality: \[ T_N = \frac{\widehat \theta_N - \theta_0}{\mbox{s.e.}(\widehat \theta_N)}, ~~~~~~~~~~~ T_N \stackrel{d}{\rightarrow} N(0,1) ~~ \mbox{under} ~ H_0\]

Test Using Asymptotic Normality

• Compute asymptotic p-value based on asymptotic normality: \[ 2 (1 - \Phi(| T_N | ))\]

• To test at level \(\alpha\), reject if p-value less than \(\alpha\).

Example

\[\ln(wage)_i = \beta_0 + \beta_1 \mbox{exp}_i + \beta_2 \mbox{exp}_i^2 + e_i\]

Example

Using studentized test statistic and asymptotic normality: \[ T_N = \frac{\widehat \beta_2}{\mbox{s.e.}(\widehat \beta_2)}, ~~~~~~~~~~~ T_N \stackrel{d}{\rightarrow} N(0,1) ~~ \mbox{under} ~ H_0\]

test.stat <- coef(reg.1)[3]/sqrt(vcovHC(reg.1)[3,3])
p.value1 <- 2 * (1- pnorm(abs(test.stat)))

[1] "Test Statistic: -2.89"

[1] "P-Value: 0.00391"

Bootstrap Hypothesis Test

Again consider

• \(H_0: ~ \theta = \theta_0,\)
• \(H_1: ~ \theta \ne \theta_0\).

Again define studentized test statistic as: \[ T_N = \frac{\widehat \theta_N - \theta_0}{\mbox{s.e.}(\widehat \theta_N)}\] computed on original sample.

Bootstrap Hypothesis Test

Again define studentized test statistic as: \[ T_N = \frac{\widehat \theta_N - \theta_0}{\mbox{s.e.}(\widehat \theta_N)}\] computed on original sample.

• Use same studentized test statistic as before, but use bootstrap distribution of studentized test statistic instead of normal distribution.

Bootstrap Hypothesis Test

• Given \(B\) bootstrap replications, compute p-value
  \[\frac{1}{B} \sum_{i=1}^B 1\{ |T^*(b) | > |T_N| \},\] where \(T_N\) is test statistic on original sample and \[T^*(b) = \frac{\widehat \theta^*(b) - \widehat \theta}{s(\widehat \theta^*(b))}.\]

• To test at level \(\alpha\), reject if p-value less than \(\alpha\).

Example

• \(H_0: ~ \beta_2 = 0,\)
• \(H_1: ~ \beta_2 \ne 0\).

Define:

\[ T_N = \frac{\widehat \beta_2 }{\mbox{s.e.}(\widehat \beta_2)}\] \[ T^*(b) = \frac{\widehat \beta_2^*(b) - \widehat \beta_2}{s(\widehat \beta_2^*(b))} \]

Example

f.boot.p <- function(d,i){ 
  reg.b <- lm(I(log(wage))~ experience+ I(experience^2),data=d[i,])
  test.stat <- (coef(reg.b)[3]-  coef(reg.1)[3])/sqrt(vcovHC(reg.b)[3,3])
  test.stat
  }
boot.test<- boot(data=df, statistic=f.boot.p, R=B, stype="i")
p.value.boot <- mean(abs(boot.test$t) >= abs(test.stat))

Example

Use studentized test statistic but bootstrap distribution of studentized test statistic.

f.boot.p <- function(d,i){ 
  reg.b <- lm(I(log(wage))~ experience+ I(experience^2),data=d[i,])
  test.stat <- (coef(reg.b)[3]-  coef(reg.1)[3])/sqrt(vcovHC(reg.b)[3,3])
  test.stat
  }
boot.test<- boot(data=df, statistic=f.boot.p, R=B, stype="i")
p.value.boot <- mean(abs(boot.test$t) >= abs(test.stat))

[1] "P-Value: 0.005"

Comparing

Paired vs Residual vs Wild Bootstrap

• We have done nonparametric bootstrap for OLS, also called “pairs bootstrap” in regression context.

• Other bootstrap procedures for regressions:

  • Residual bootstrap,

  • Wild bootstrap.

• Tradeoffs?

Non-i.i.d data

• We have done bootstrap for i.i.d. data.

• Can modify bootstrap procedure for

  • clustered sampling
    • straight forward modification, resample clusters
  • time series
    • less straight forward (block bootstrap).